#include <stdlib.h>
#include <stdio.h>
#include <math.h>
/* take as argument the base filename and output .avg file with averages of the first two columns; the third column will the the standard deviation of x and the fourth will be sd of y
 */
int main(int argc, char *argv[]) {
  if (argc<3||argc>4) {
    printf("Usage: avg filebase nfiles\n");
    return;
  }
  else {
    FILE *fp;
    char fname[150], string[500];
    float col1[200][200], col2[200][200]; // expect no more than 100 data points and no more than 100 files
    float sum1[200], sum2[200],var1[200], var2[200], avg1[200], avg2[200], sum=0, summ=0;
    int i,j, nfiles, ndata=0;
    int ndatamin=1000;
    nfiles = atoi(argv[2]);
    
    for (i=0;i<200;i++){sum1[i]=0;sum2[i]=0;var1[i]=0;var2[i]=0;avg1[i]=0;avg2[i]=0;}
    for (i=0; i<nfiles; i++) {
      // open each file and start computing the sum
      sprintf(fname, "%s.%d",argv[1],i); //check this later
      if (argc==4) sprintf(fname, "%s.%d.%s", argv[1], i, argv[3]);
      //printf("filename: %s\n", fname);
      fp = fopen(fname, "r");
      ndata=0;
      while (!feof(fp)) {
	if (fgets(string, 499, fp)){
	  sscanf(string, "%g %g %*g", &col1[i][ndata], &col2[i][ndata]);
	  sum1[ndata]+=col1[i][ndata];
	  sum2[ndata]+=col2[i][ndata];
	  ndata++;
	}
      }
      if (ndata<ndatamin) ndatamin=ndata;
      fclose(fp);
    }
    /* now compute the average for each row */
    for (j=0; j<ndatamin; j++) {
      avg1[j]=sum1[j]/nfiles;
      avg2[j]=sum2[j]/nfiles;
    }
    /* now compute the variance for each row */
    for (j=0; j<ndatamin; j++) {
      sum=0; summ=0;
      for (i=0; i<nfiles; i++) {
	sum+=(col2[i][j]-avg2[j])*(col2[i][j]-avg2[j]);
	summ+=(col1[i][j]-avg1[j])*(col1[i][j]-avg1[j]);
      }
      var2[j]=sqrt(sum/nfiles);
      var1[j]=sqrt(summ/nfiles);
    }
    
    /* write the output */
    printf("ndata:%d\n", ndatamin);
    sprintf(fname, "%s.avg",argv[1]);
    fp = fopen(fname, "w");
    for (j=0; j<ndatamin; j++) {
      fprintf(fp, "%g %g %g %g\n", avg1[j], avg2[j], var1[j], var2[j]);
    }
  }
}
